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Fourier methods for the perturbed harmonic oscillator in linear and nonlinear 

Schrodinger equations 

Philipp Badeo and Sergio BlaneqJ 
Universitat Politecnica de Valencia, Instituto de Matemdtica Multidisciplmar, E-46022 Valencia, Spain 

We consider the numerical integration of the Gross-Pitaevskii equation with a potential trap given 
by a time-dependent harmonic potential or a small perturbation thereof. Splitting methods are 
frequently used with Fourier techniques since the system can be split into the kinetic and remaining 
part, and each part can be solved efficiently using Fast Fourier Transforms. To split the system into 
the quantum harmonic oscillator problem and the remaining part allows to get higher accuracies in 
many cases, but it requires to change between Hermite basis functions and the coordinate space, 
and this is not efficient for time-dependent frequencies or strong nonlinearities. We show how to 
build new methods which combine the advantages of using Fourier methods while solving the time- 
dependent harmonic oscillator exactly (or with a high accuracy by using a Magnus integrator and 
an appropriate decomposition). 

PACS numbers: 02.60.-x, 02.60.Cb, 02.60.Lj, 02.70.Hm 
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I. INTRODUCTION 

The numerical integration of the Gross-Pitaevskii 
equation (GPE) 



9 ,, s 



i 

2/1 



A + V(x, t) + cr(t)\i/)(x, t)\ 2 ij)(x, t) 



x £ M. d , describing the ground state of interacting bosons 
at zero temperature, the Bose-Einstein condensates, has 
attracted great interest [lrl3] after the first experimental 
realizations [4| . We present a new efficient way to solve a 
special class of GPE, namely that of weakly interacting 
bosons in a single time-dependent trap. To be more spe- 
cific, the potential trap V is taken to be a perturbation of 
the (time-dependent) d-dimensional harmonic oscillator, 
i.e. V(x,t) = x T M{t)x + eV!{x,t) where M(t) e R dxd 
is a positive definite matrix and eV~i(x,t) is a small per- 
turbation. The real scalar function a originates from the 
mean-field interaction between the particles and corre- 
sponds to repulsive or attractive forces for positive or 
negative values of a(t), respectively J5J. Notice that the 
non-interacting case, a = 0, corresponds to the linear 
Schrodinger equation. 

Several methods have been analyzed to compute both 
the time evolution and the ground state of the GPE in 
the course of the last decade [IH3|,[6|,0]> among them finite 
differences, Galerkin spectral methods and pseudospec- 
tral methods for Fourier or Hermite basis expansions. It 
has been concluded Q that these pseudospectral meth- 
ods perform best for a wide parameter range for the GPE. 
The Fourier type methods can be implemented with Fast 
Fourier Transform (FFT) algorithms since the trapping 
V causes the wave function to vanish asymptotically thus 
allowing to consider the problem as periodic on a suffi- 
ciently large spatial interval. Their advantages are high 



accuracy with a moderate number of mesh points and 
low computational cost. For harmonic oscillator (HO) 
problems, however, the exact solution is known and by 
expanding the solution in Hermite polynomials, highly 
accurate results are obtained if the HO is solved sepa- 

ratelyiliLlL 

It is claimed [2|, that either Hermite or Fourier pseu- 
dospectral methods are the most efficient, the choice de- 
pending on the particular parameter set. Motivated by 
these results, we show how both methods are combined to 
retain both the accuracy of the Hermite method and the 
speed of the Fourier transforms, i.e. to rewrite the Her- 
mite method as a single simple pseudospectral Fourier 
scheme. We have found, that this approximation per- 
forms, for the studied problem class, always equal to or 
better than the original Fourier method and therefore 
has to compete with Hermite expansions only. Hermite 
schemes suffer from large computational costs when the 
number of basis terms in the expansions is altered along 
the integration or taken very large, being the case for 
time dependent trap frequencies M(t) or strong nonlin- 
earities a(t). It is in this setup, where our new method 
substantially improves the Hermite performance, and it 
can indeed be regarded as the optimal choice for the num- 
ber of Hermite basis functions (for an equidistant grid) 
at each time step. 

For the ease of notation, we restrict ourselves to the 
one-dimensional problem 



ijgl> = H (t)1> + {eVi(x, t) + <T(t)\Tjj\ 2 ) 



1> 



where 



Ho(t) = ^-p 2 + \(jJ 2 {t)x 2 



(1.1) 



(1.2) 
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and p = —i-^- The boundary conditions imposed by the 
trap require the wave function to go to zero at infinity, 
and up to any desired accuracy, we can assume ip(x,t) 
and all its derivatives to vanish outside a finite region, say 



[a, b], which we divide using a mesh (usually with N — 
2 k points to allow a simple use of the FFT algorithms). 
Then, the partial differential equation (jl.lj) transforms 
into a system of ordinary differential equations (ODEs) 

*—«(*) - H (t)u(t) + (eVi(X,t) + a(t)\u(t)\ 2 ) u(t), 

(1.3) 



and the harmonic part becomes 



1 



H (t)=T + ^ 2 (t)X 2 , 



(1.4) 



with u € C*^, where Uj(t) ~ ip(xi,t), Xi — a + ih, i = 
0, 1, . . . , N — 1, h = (b — a)/N, X = diag{xo, . . . , xjv-i} 
and T denotes a discretization of the kinetic part. 

This system of nonlinear ODEs can be numerically 
solved by standard all purpose ODE methods. However, 
because of the particular structure of this problem, dif- 
ferent numerical methods can differ considerably in ac- 
curacy as well as in computational cost and stability. In 
addition, the structural properties of the system lead to 
the existence of several preserved quantities like the norm 
and energy (for the autonomous case). 

The accurate preservation of these quantities as well as 
the error propagation and performance of splitting meth- 
ods explain why they are frequently recommended for the 
time integration [2|, y, l6| and make them subject of in- 
vestigation in this work. 



II. SPLITTING METHODS 

Let us consider the separable system of ODEs 

v! = A(u)+B(u), u(t ) = uo (= C N , (11.1) 

where we assume that both systems 

v! = A(u), v! = B(u) (11.2) 

can either be solved in closed form or accurately inte- 
grated. If ip\ , tp\ represent the exact flows associated 
to (III. 21) then, to advance the solution one time step, h, 
we can use, for example, the composition ijj h = ip h oip^ 

which is 



(vjfW)), 



(i.e. u(t + h) ~ ^M = ¥>L A1 

known as the first-order Lie- Trotter method. A method 
has order p if ip^' = tph + 0(h p+1 ) where ipt denotes 
the exact global flow of (|II.1[) . Sequential application of 
the two first-order methods Vv an< ^ ^ ts adjoint ?/v ~ 
<p\ 'o^' with half time step yields the second order 
time-symmetric methods 



,[2] [A] [B] [A] 

Vh,A = <Ph/2 °<Ph ° Vh/i 



Vh,B 



[B] [A] [B] 

<Ph/2°^h °Vh/2 



(11.3) 

(II.4) 



(referred as ABA and BAB compositions) . The contrac- 
tion via the (l-parameter-)group property of the flows 



that eliminated one computation is called First Same As 
Last (FSAL) property and can also be used with higher 
order m— stage compositions 



/ [p] W [B] 



[A] [B] 

, ' 0( Paih 0( Pb 1 h 



(11.5) 



if a m = or b m = and repeated application of the 
scheme without requiring output. For linear problems, 
it is usual to replace the flow-maps by exponentials (for 
nonlinear problems the same is possible using exponen- 
tials of Lie operators). In this notation, the equation 
iu' = A(u) + B(u), whose formal solution for the evolu- 

is ap- 



tion operator is denoted by (jr t = e - lt ( A + B ) . 

proximated for one time step, h, by the order p compo- 
sition pi.5p or, equivalently, 



„/,[p] — „-iha m A -ihb m B 



— iha\A —ihb±B 



(II.6) 



We keep in mind that, in a nonlinear problem, if B 
depends on u, it has to be updated at each stage because 
u changes during the evolution of er lha i A . 

Backward error analysis shows that the action of a 
splitting method is equivalent to solve exactly, for one 
time step h, a perturbed differential equation (for suffi- 
ciently small h) that can contain higher order commuta- 
tors 



= (A + B + hai t i[A,B] 

+h 2 (a 2A [A,[A,B]}+a 2 ,2[BdB,A]}) 



(II.7) 



([A, B] :— AB — BA, and for simplicity, we have denoted 
A = A(u) ■ V, B = B(u) ■ V) where the coefficients a.;, bi 
are chosen in order to cancel out the coefficients a^j up 
to a given order. It is thus important to analyze the dom- 
inant contributions to the error from the commutators, 
in order to choose the most appropriate method or to 
build new ones. 

There exist many different splitting methods which are 
designed for different purposes, depending on the struc- 
ture of the problem, the desired order, the required sta- 
bility, etc. g-tLJ]. If the operator A corresponds to the 
kinetic part, which is quadratic in momenta, and the op- 
erator B is the potential, diagonal in coordinate space, 
the commutator [B, [B, [B,A]]] vanishes and partitioned 
Runge-Kutta-Nystrom methods become favorable. The 
coefficients for this family of methods have to solve a 
significantly reduced number of order conditions and, in 
general, their performance is superior to splitting meth- 
ods designed for general separable problems, in addition 
to being more stable. Efficient schemes of order 4 and 6 
are obtained in (^]. On the other hand, when Hq is the 
dominant part, it is worth to take a closer look at the 
split (jl.ll) . This would correspond to ||5|| <C ||A|| and for 
this case, when facing autonomous problems, there exist 
tailored methods which have shown a high performance 
in practice. Writing the equation as iu' = (A + eB)u 
(with e a small parameter), it is clear that the local er- 
ror of the second order methods (|II.3|) or (|II.4[) comes 
from the commutators at third order ([A, [A, eB]] and 



[eB, [A, &B]J) and we can say that the local error is of 
order 0(eh 3 + e 2 ft, 3 ). The coefficients Oj, 6j in the general 
composition (|II.6|I can be chosen to cancel the dominant 
error terms, say, the 0(eh r ) terms for relatively large val- 
ues of r. Then, one can denote the effective order of a 
method by (r, p) with r > p when the local error is given 
by 0(eh r+1 + e 2 hP +1 ). The method is of order p, but 
in the limit e — ► it is considered as of order r > p. 
Using this split allows to gain a factor e in the accu racy 
even for general splitting methods where r = p. In [111 ] , 
several methods of order (r, 2) for r < 10 are obtained 
with all coefficients a% , bi being positive and some other 
schemes of order (r, 4) for r = 6, 8 are presented. For 
near-integrable systems, these last methods are the most 
efficient and stable ones. Despite the gain of accuracy, 
the split into a dominant part and a small perturbation 
is left unconsidered when it leads to involved or compu- 
tationally costly algorithms. This issue is addressed in 
this work. 

To take account for the time dependence in the vector 
fields in (11.11) . a more detailed analysis is required. The 
general separable equation 



H of |LT|). 



iu' = A(u,t) + eB(u,t) 



(11.8) 



can be solved by considering the time as two new inde- 
pendent coordinates 



iu' = A(u,t2) 
t'i = 1 



= e£(u,ti) 



t' = 1 



(II.9) 



d I 



^-p 2 + ^ficj 2 {t)x 2 ) r 



(11.11) 



is easily computed for a time step using Fourier trans- 
forms. Before giving the details on the time integration, 
some remarks on the formal solution are necessary. 

It is well known that Hq is an element of the Lie alge- 
bra spanned by the operators {E = x 2 /2, F — p 2 /2, G = 
\ (px + xp) } , where \i = 1 for simplicity, and its commu- 
tators are 

[E, F] = iG, [E, G] = 2iE, [F, G] = -2iF. 

This is a three-dimensional Lie algebra and the solution, 
ip(x,t) = U(t,0)ip(x,0), of pi.ll[) can be expressed as a 
single exponential using the Magnus series expansion [16J, 
LLZ| or as a product of exponentials [18| ■ It is possible to 
formulate the evolution operator U(t, 0) in many different 
ways, the most appropriate depending on the particular 
purpose, e.g. using the Magnus expansion 

U(t,0) = exp{f 1 (t)E + f a (t)F + f 3 (t)G) (11.12) 



for certain functions /»(£) [19| . Approximations of (|IL12|) 
for one time step, h, on the other hand, are easily ob- 
tained, e. g. a fourth-order commutator-free method is 
given by |20J 



This standard split is equivalent to a new system in an 
extended phase space which is not near-integrable and 
the highly efficient and stable methods designed for those 
problems lose their excellent performance. Following [15j . 
the near integrability is recovered if we introduce only one 
variable as follows 



in 



A(u,h) 
l 



= eB(u,h). (11.10) 



This split requires to solve exactly the equation iu' = 
A{u, t) for the corresponding fractional time steps and to 
freeze the time when solving iu' = eB(u, t±) (see Tabled]). 



TABLE I. Algorithm for the numerical integration of the sys- 
tem (|II. lOf) by the composition (|H.5|) . 



te[t li - 1] ,t li - 1] + bih] 

te[t [i - 1] ,t^- 1] +aih] 



Algorithm for one time step t n — > t n + h 

z — t n 
do i — 1, m 

solve : u = B(u,t [i ~ 1] ), 
solve : u' = A(u,t), 

t [i] =t [t - 1] +a t h 
enddo 

+ f[m] 

f-n + 1 — v 

We show that the exact solution of the non- 
autonomous problem, in our setting the dominant part 



U(t + h,t) 



LI 

2 ( 2 J 

h,\ 



ex Pl o(o? + w Lo X ) 



(11.13) 



.2 1 „2- 



xexp[--(-p 2 +^ R -x'))+0(h») 



where 



2 2 i a 2 

uj l = au; 1 + puj 2 , 



= Pu 2 



= uj(t n + ah), d = | - &, c 2 



2^6' 

can be considered 



with u. 

and a ~ ■k — -y=, 13 — 1 — a. It 

as the composition of the evolution for half time step 

of two oscillators with averaged frequencies, using the 

fourth-order Gauss-Legendre quadrature rule to evaluate 

u>(t). Different quadrature rules can also be be used and 

correspond to different averages along the time step, see 

[13,|20|. In the limit when u is constant the exact solution 

is recovered. Higher order approximations are available, 

if more accurate results are desired, by approximating the 

functions fa in (III.12[) via truncated Magnus expansions. 

Our objective is to obtain a factorization of the solu- 
tion which only involves terms proportional to E or F 
since they are easy to compute as we show in the follow- 
ing paragraph. 

Starting from pi.!3[) or high order approximates of 
pi.!2|) , the main result of this work is the natural decom- 
position for the application of Fourier spectral methods. 



A. Hamiltonians for spectral methods 

We now analyze how to compute the evolution of dif- 
ferent parts of the Hamiltonian by spectral methods. 
The spatial derivative (or kinetic part) associated to the 
scmidiscretized problem (jOJ can be solved in the mo- 
mentum space by noting that 



i— u{t) = Tu{t) = F N D N T N u, 



(11.14) 



where Dn is diagonal and Fn denotes the discrete 
Fourier transform of length N, whose computation can 
be accomplished by the FFT algorithm with 0(N log N) 
floating point operations. The solution of (|II.14[) for one 
time step, h, is given by 

u(t + h) =T N 1 e- lhL>N T N u(t) 

which requires two FFT calls. The exponentials in 
e " need to be computed only once and can be reused 
at each step such that the cost of the action of e~ lhDN 
corresponds to N complex products. 

For the remaining part, the following well-known result 
is very useful 

Lemma II. 1 If F is real valued, the equation 

i—<f>(x, t) = F(x, \<t>(x, t)\)4>(x, t), (11.15) 

leaves the norm invariant, \<f>(x,t)\ = \<$>(x, 0)|, and then 

tj>(x, t) = e- itFix '^ xm (f>(x, 0). (11.16) 

On the other hand, it is well known that the solutions 
of the linear Schrodinger equation with the harmonic po- 
tential 

i^{x,t) = ^{p 2 + x 2 )(j>{x,t) (11.17) 

can be expressed in terms of Hermite polynomials: 

<ftx,t) = ^c n e- iE - t h n {x) (11.18) 

n=0 

where 

E n =n+ -, hJx) = — - — ■ H n (x)er x2/2 

(11.19) 
and H n (x) are the Hermite polynomials satisfying the 
recursion 

H k+1 (x) = 2xH k {x) - 2kH k _ x , k = 1, 2, . . . 

with Hq(x) = 1, Hi{x) = 2x. The weights c n 
can be computed from the initial conditions, c„ = 
J h n (x)4>(x,0) dx. 



The previous results show how to compute individual 
parts of the equation and thus permit different ways of 
splitting the system in two solvable parts 



iifH = (A + B)tp 

We consider the following cases: 
(i) Fourier(F)-split. 



(11.20) 



A 



1 



;i r 



B(t) 



u(t) 



x 2 +eV I (x,t)+a(t)\iP\ 2 , 



2" ' w 2 

(11.21) 

We take the time as a new coordinate, as in (jll.101) . and 
evolve it with A, which is now autonomous and exactly 
solvable, and freeze the time in B, which is then solved 
using the result from Lemma III. 11 Here, A and B are 
diagonal in the momentum and coordinate spaces, re- 
spectively, and we can change between them using the 
Fourier Transforms. 

(ii) Harmonic oscillator(HO)-split. Let for a mo- 
ment w = 1 , the Hermite expansion then suggests a split 

A = ^(p 2 +x 2 ), B(t)=eV I (x,t) + a(tM 2 , (11.22) 

where the solution for the equation ivl — Au can be ap- 
proximated using a finite number of Hermite basis func- 
tions, i.e. 



<f> M (%,t) 



M-\ 

E 

71 = 



( ,, t 



-iE„t 



h n (x). 



flL22p) 



Since B is diagonal in the coordinate space it will act 
as a simple multiplication if we choose a number of Her- 
mite basis functions and evaluate them at the points of 
a chosen mesh (e.g. using the Gauss-Hcrmite quadrature 
[21j or on equidistant grid points [2j). This split can be 
of interest if all contributions from B(t) are small with 
respect to A and methods for near-integrable systems are 
used. If the frequency is time dependent, the correspond- 
ing split is 



A{t)= 1 -{p 2 +u 2 {t)x 2 ), 



B{t) 



eV I (x,t)+a(tM 2 . 

(11.23) 
In this case, it is convenient to take the time, t as a 
new coordinate as shown in (III.lOp . The solution of 
ivl = A(t)u, in the algorithm in Table HI can be approx- 
imated by the Magnus expansion, e.g. (III.12[) or (III. 131) . 
but these factorizations are not appropriate for use with 
spectral methods since it would require two sets of basis 
functions and additional transformations. 

In general, the split (i) can be considered faster and 
simpler since A = T can be computed in the momentum 
space, and one can easily and efficiently change from mo- 
mentum to coordinate space via FFTs. The choice (ii), 
on the other hand, allows us to take advantage of the 
structure of a near-integrable system if, roughly speak- 
ing, ||S|| < ||A||, but it requires to solve the equation 



for the (time-dependent) harmonic potential exactly (or 
with high accuracy) . The evolution of the constant oscil- 
lator is easily computed using Hermite polynomials (see 
@,Sl2it)- As mentioned, methods tailored for this family 
of problems show a high performance. The main draw- 
back, however, is that the evolution for e ~ tha i A has to be 
carried out using a basis of Hermite polynomials whereas 
the B part is advanced in space coordinates 2l| resulting 
in possibly costly basis transformations. 



B. Solving the Harmonic oscillator by Fourier 
methods 



We propose a new method which combines the advan- 
tages of both splittings. It retains the advantages of the 
HO-split (ii) while being as fast to compute as the F- 
split in (i). For this purpose, we briefly review some 
basic concepts of Lie algebras. 

Given X, Y two elements of a given Lie algebra it is 
well known that 



„x v -x 



1, 



Ye~* = Y + [X, Y] + -[X, [X, Y}] + ... 

If {Xi, . . . , Xk} are the generators of a finite dimensional 
Lie algebra and 



k 

E 

71=1 



oi n X n , G 



k 

E 

71=1 



Qn^n 



then 



G(i)=e tz Ge- tz = ^a„(i)X n 

n=l 



(11.24) 



where G(t) satisfies the differential equation 

G'(t) = [Z,G(t)}, G(0) = G. 

Notice that if G = x or G — p and Z = H(x,p) these 
equations are equivalent to the equation of motion for the 
classical Hamiltonian H(x,p). It is also well-known that 
two operators are identical on a sufficiently small time 
interval if they satisfy the same first order differential 
equation with the same initial conditions 22] . 

Given X(x),P(p),W — \{p 2 + x 2 ) and an analytic 
function F(x,p), we are interested in the following ad- 
joint actions 



e- ltx F{x,p)e ltx =F(x,p + tX') 
e- ltp F{x,p)e ltp = F(x - tP',p) 
e- ltw F{x,p)e ltw = F(x,p) 

where X' = dX/dx, P' = dP/dp and 

x = cos(t)x + sin(t)p 
p = — sin(t)x + cos(t)p. 



(II.25a) 
(II.25b) 
(II.25c) 



In classical mechanics, this corresponds to a kick, a drift 
and the harmonic oscillator rotation. 

As we have seen, the exponentials can be easily com- 
puted by Fourier spectral methods. It is then natural to 
ask the question if it is possible to write the solution of 
(|II.11|) as a product of exponentials which are solvable 
by spectral methods. The answer is positive and it is 
formulated as follows. 

Let us first consider the pure harmonic oscillator, 
whose result was obtained in [23j], and for which we 
present a new proof that applies equally to the general 
case. 

Lemma II. 2 Let A\ = ^p 2 , B\ = ^x 2 and 

g(t)=ma(t), f(t) = tan(t/2) . (11.26) 

Then, the following property is satisfied for \t\ < n: 



-it(Ai+Bi) _ e -if(t)Ai e -ig{t)Bi & -if{t)Ai 

- P ~if(t)Bi p -ig{t)A x „-if(t)Bi 



(11.27) 

(11.28) 



Proof. It suffices to impose that the right hand side of 
(|TL27)) (or (|IL"2"8"J) ) satisfies the equation pL17l) with the 
identity as initial condition what can be verified by differ- 
entiating and simple applications of the rules pi.25al b). 
A more constructive way to derive the functions /, g 
makes use of the parallelism with the one-dimensional 
classical harmonic oscillator with Hamiltonian H = 



b 2 



1„2 



q and Hamilton equations 




= {A 




where 



.4 



1 




B 




(11.29) 



(11.30) 



The Lie algebra generated by the matrices A, B is the 
same as the Lie algebra associated to the operators 
A\ , B\ for the Schrodinger equation with the harmonic 
potential pi. 170 . 

The exact evolution operator of pl.29j) is 



0(t) 



cos(t) sin(t) 
— sin(t) cos(t) 



(11.31) 



which is an orthogonal and symplectic 2x2 matrix. For 
the splitted parts, the solutions are easily computed to 



J(t)A _ / 1 /(*) 



pfl(*)-B _ 



1 







\o i ; \-g(t) i 

and then, equating the symmetric composition 



e fA e gB e fA = l 1-f-g 2f-f-g 



l-f-9 



to (|II.3ip . we obtain (III.26|) which is valid for \h\ < n. 
The decomposition (|II.28|> is derived analogously. Using 
the Baker-Campbell-Haussdorf-formula, it is clear, that 
both results remain valid, up to the first singularity at 
t = ±7r, when replacing the matrices A, B by the corre- 
sponding linear operators A\ , B\ , since all computations 
are done in identical Lie algebras. □ 



Then, the following decomposition 



e -iI(i„ s +U iix 2 ) e -if{i/+^l, 3 ) 
e -if(t)ix 2 e -ig(t)±p 2 e -ie(t)±x 2 



(11.36) 



is satisfied for < t < t* , where t* is the smallest positive 
root of g(t). 



It is immediate to generalize this result to the equation 

,2 



i—ip(x,t) 



p 2 + /iyX 2 j r(.rJ). 



/j, u) > by replacing (|II.26|) with 

g= — sin(wi), f = fiuj t&n ( —t) . 

[ibJ V 2 / 

This result is valid for \t\ < t* = ir/ui. 

The following theorems extend this idea to decompo- 
sitions of operators that appear after the approximation 
of the time dependent parts via (|II.12|) or by the compo- 
sition (|rLl3|) . 



Theorem II. 3 Let a, /3,7 be constants, r\ = ^Jorf—JP 
and 

g(t)= 1 /r ] -sm(r 1 t), (11.32) 

•f W = -J7\( 1 - cos ^ + ~ sin ^) 

g(t) V v 

e(<) = -77s ( 1 - cos(r/t) - - sm(nt) 
Then, the following decomposition holds for < t < n/n: 

— i^(ax 2 +fl(xp+px)+^p 2 ) 

= e -if(t)^x 2 e ~ig{t)hp 2 e ~ ie W2 x2 . (11.33) 

Proof. The proof follows the lines of the proof of 
Lemma III. 21 The evolution operator associated to the 
classical Hamiltonian H = ^ (ax 2 + 2/3xp + 7p 2 ) is given 
by 

/ cos(r)i) + & sin(r?t) 2 sin(?7£) 

V -% sin(ryi) cosfai) - | sin(?yt) _ 

and equality to the right hand side of (|II.33j) is verified by 
straight-forward computation of the matrix exponentials. 
The solution is valid until the first singularity at t = ir/n. 
Using pi.25a[ b). it can be checked that both sides of 
(|II.33| satisfy the same differential equation and initial 
conditions. □ 

Theorem II. 4 Let ujk £ R, Ck — cos(ojfet/2), s/, — 
sin(o;/ £ t/2) for k = L, R and 



(11.34) 



g(t) = SlCr/ul + clSr/ur, 
/(*) = ~77V ( l ~ ClCr + — s lSr 

g(t) V u R 

I \ 1 ( Ur 

e(t) = —r- 1 - c L c R H slsr 

g(t) V u L 



(11.35) 



The proof is similar to the previous one. 



III. THE HERMITE-FOURIER METHODS 

With the presented exact decompositions at hand, we 
now solve the discretized GPE (|I.3|) by splitting methods 
using the symmetric compositions (111.31) and (|II.4|) . Let 
us first consider the case w = 1 and take the HO split 
A = A 1 + B 1 



^ [2] 
Vh,B 



e -ih{A 1 +B 1 )/2 e -ihB e -ih(A 1 +B 1 )/2 (ffl.l) 

e -ihB/2 e -ih(Ai+Bi) e -ihB/2 



{111.2) 



Repla cing the exponentials e ~ ih( - Al+B ^ by pi.27|) or 
pi.28|) , we obtain four different methods whose computa- 
tional costs differ considerably. At first sight, using the 
FSAL property, both (|III.1[) and (IIII.2I) are equivalent 
from the computational point of view and require one 
exponential of B and another one of A\ + B\ per step. 
However, a significant difference arises when we plug in 
the decompositions (|II.27|) or (jll.28[) . Only the combi- 
nation (JIII.2J) with pi.28|) yields a method that involves 
only a single FFT call per step 



< 



[2] _ e -ihB/2 -ifc(Ai+Bi) e -ihB/2 

-ihB/2 _-i/(h)Bi e -ig(h)Ai „-i/(fc)Bi „-ihB/2 



e -i(hB/2+f(h)Bi) e -ig(h)Ai e -i(hB/2+f(h)Bi) 



(111.3) 

and solves exactly the harmonic oscillator for \h\ < h* . 
For any other combination, more kinetic terms have to 
be computed per step since the FSAL property cannot be 
exploited to full extent and hence result in more costly 
24] methods for the same accuracy. 

The general composition (|II.6[) can be rewritten in the 
same way by replacing each flow e~ laiA in (|II.6I) by the 
composition (jll.28[) 

a ^ e — i{hb m B+a m Bi) e -ig(a m h)A 1 g — i(hb m B+a m - 1 B 1 ) 

-i{hb 1 B+u 1 B 1 ) -ig(a 1 h)A 1 -ia a Bi (tjt a\ 

where a^ — f(ak+\h) + f(akh), k = 0, 1, . . . , m + 1 with 
«o = o-m+i = 0. This method is valid for \aih\ < h* , i — 
1 , . . . , m and requires only m calls of the FFT and its 
inverse, but reaches the same accuracy as if the Hermite 
functions were used. 



In the more general case with a time-dependent fre- 
quency, u>(t), starting from the HO splitting, the time- 
dependent part is first approximated by Magnus expan- 
sions pTT2|) or pTT3|) . Theorems HOI and HOI then pro- 
vide decompositions to write the product of exponentials 
in a similar way as in (|III.4j) but now au = f(ak+ih) + 
e(akh) and it is valid for \a,ih\ < t* , i = 1, . . . , m where 
t* is the first zero of g(t). At each stage, one has to 
compute f(a,ih),g(aih),e(aih) from ui{t). 

As in the previous case, it only requires m calls of 
the FFT and its inverse, like the standard Fourier pseu- 
dospectral methods. For stability reasons, it seems 
convenient to look for splitting method whose value of 
maxi{|a.;|} is as small as possible. 




-7T -2 



2 7T 3 



FIG. 1. Error in logarithmic scale for the integration of the 
ground state of the Harmonic potential using the split (|H.28|) 
for T £ [— 7r,7r] (integration forward and backward in time). 
The left and right panels show the 2-norm error and a zoom 
about T — n, respectively. 



IV. NUMERICAL EXAMPLES 

We analyze first the performance of the methods con- 
sidered in this work for the one-dimensional problem (|I.1|) 
with (jl = w = 1, and the pure harmonic trap, i.e. 
Vi = 0. 

To illustrate the validity of the decomposition pre- 
sented in Lemma III. 21 we first consider the linear prob- 
lem (a = 0). We take the ground state at t = as initial 
condition whose exact solution is given by 



ip(x,t) = 



1 



7r 1 /4 % /2"n! 



-.e^^e- 1 



72 



We discretize on the interval [—10, 10], to ensure the wave 
function and its first derivatives vanish up to round off 
at the boundaries, and sample it at N = 1024 equidis- 
tant grid points. We integrate, with only one time step, 
from t = to T for T e [— 7r,7r], i.e. forward and back- 
ward in time. We measure the integrated error in the 
wave function, ||w ea; (r) — w Q (T)||2, where u a (T) denotes 
the approximate numerical solution obtained using the 
split (JII.28J) and u ex (T) is the exact solution at the dis- 
cretized mesh. The result of this comparison is illustrated 
in Fig. Q] (left). The split pL28t reproduces, for \T\ < n, 
the exact solution up to round off, as expected. The 
right panel in Fig. [1] displays a zoom near a singularity 
where the error grows rapidly due to double precision 
arithmetic. 

We analyze how the approximation properties of the 
Hermite decomposition (|II.22b|) strongly depend on the 
function in question and on the chosen number of ba- 
sis functions, M. We compute the M required to 
reach round-off precision for the evolution of a dis- 
placed ground state as initial condition, ips(x,0) = 
e -(z-<5)72/( 7r i/4 v /2^i) f rom t = to T = 10 in one 
time step. From initial conditions computed on a mesh, 
this can be accomplished as follows [2] : 



Uex(T) 



-iT{A l +B 1 ), 



u - K 1 e- li Ul K u Q (IV. 1) 



where Z?i 



diag{ 



1 3 



2M-1" 



., . ,, ., -}, M is the number of 

basis elements considered, and K^j — hi-i(xj), i = 



1, . . . ,N = 512 with h n (x) given in (|jLT8l) . 

For 5 - 1 



l,...,M,j 

x G [—10, 10]. For 8 = j^j, round off accuracy is achieved 
with M = 8 while for i5 = 2 it is necessary to take 
AI = 29. We observe that the Hermite decomposition 
is very sensitive to the initial conditions [21] . The Her- 
mite basis works efficiently as far as the initial conditions 
as well as the evolution thereof can be accurately approx- 
imated using a few number of basis elements and one has 
to keep in mind that, for nonlinear problems, the number 
of basis functions necessary to reach a given accuracy can 
vary along the time integration. 

Next, we study the following values for the nonlinear- 
ity parameter: a = 10~ 2 ,1,10 2 . The case a = 10~ 2 
illustrates the performance of the new methods if ap- 
plied to problems (linear or nonlinear) which are small 
perturbations of the Harmonic potential whereas the val- 
ues (7 = 1, 10 2 are large enough to demonstrate the non- 
linearity effects on the approximation properties of the 
methods. Physically, a is proportional to the number of 
particles in a Bose-Einstein condensate and to the inter- 
action strength [5]. 

For all cases, we choose the initial condition ifj(x,0) — 
pe - ^" 1 ) / 2 , with p a normalizing constant. We show in 
Fig.[5]the value of |^>(a;, i)| 2 at the initial and final time. 
The spatial interval is adjusted to ensure the wave func- 
tion vanishes (up to round off) at the boundaries, here 
[-20, 20] for a = 0.01 and [-30, 30] for both a = 1 and 
a = 100 where the wave function moves faster (we only 
show the interval x G [—5,5]). One can appreciate that 
for strong nonlinearities the wave function can consid- 
erably penetrate the potential barrier and one expects 
that an accurate approximation of these wave functions 
requires a large number of Hermite functions when using 
PV.1I) and hence renders this procedure inappropriate. 



A. HO-split versus F-split 

We analyze now the advantages of the HO-split versus 
the F-split as given in (|II.23|) and (|II.21|) . 

In a first experiment, we fix the symmetric second or- 
der BAB composition (|III.2[) and apply it for both splits. 
For the HO-split, we compute the harmonic part either 
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FIG. 2. Exact evolution at t — T (solid line) from the initial 
conditions given by ip(x,0) — pe - ^ -1 ^ ' 2 (dashed line). The 
number of grid points is given by N . 



with the decomposition (|II.28|) or in the Hermite basis 
pV.ip with different numbers of basis terms. 

The parameters, initial conditions and final times are 
taken as previously for Fig. [5] Given a splitting method 
X, we denote by Xp, Xh and XhM its implementations 
with the F-split, the HO-split using the Hermite-Fourier 
method and the HO-split using M Hermite basis func- 
tions in (|IV.1I) . respectively. We measure the error ver- 
sus the number of exponentials which can be considered 
proportional to the computational cost and plot the re- 
sults in Fig. El As expected, HO-splits are advantagous 
if the system is close to a harmonic oscillator, i.e. for 
a = 0.01, 1, and if the initial conditions are accurately 
approximated by a few terms in the Hermite expansion. 
On the other hand, for strong nonlinearities a = 100, the 
Hermite polynomial based HO split shows a very lim- 
ited performance, c.f. the large number of basis terms 
in the right panel of Fig. [3] We stress that if this tech- 
nique is used for nonlinearities the number of basis terms 
should be increased along the time integration and fix- 
ing it bounds the maximally achievable accuracy and its 
limit depends on the initial condition and the strength of 
the nonlincarity. 

The Hermite-Fourier method proposed in this work 
(using the composition pi.28[l ) is clearly superior for 
weak perturbations and it keeps similar performance to 
the F-split for strong nonlinearities. 

Finally, we analyze the performance of different higher 
order splitting methods which are useful when high ac- 
curacies are desired. The following methods (whose co- 
efficients are collected in Table |TT] for the convenience of 
the reader) are considered: 

- RKNq4 (the 6-stage fourth-order method from Q). 
This is a Partitioned Runge-Kutta-Nystrom method and 
it is designed for the case where [B, [B, [B, A]}] = 0, be- 
ing the case for both the F-split and the HO-split. 

- Nh(8,2) (the 4-stage (8,2) BAB method from [TO])- 
This method is addressed to perturbed systems. One ex- 
pects a high performance if the contribution from B is 
small. 

- A/7 5 (8,4) (the 5-stage (8,4) BAB method from [HI). 
We analyze in Figures IJHH] the three problems speci- 
fied in Fig. [3J In the upper panels, the Leap-Frog meth- 



TABLE II. Coefficients for several splitting methods. 
The 6-stage 4i/i-order: SRKN 6 4 



fei = 0.0829844064174052 
b 2 = 0.396309801498368 

63 = -0.0390563049223486 

64 = 1 - 2(&i + b 2 + 63) 

65 = 63, b 6 - 62, 67 = 61 



ai = 0.245298957184271 
a 2 = 0.604872665711080 

a-i = 1/2 — (ai +02) 

04 = 03, a 5 — a.2, as = a-i 



The 4-stage (8,2) method: NI4 



61 = 1/20 

62 = 49/18 

63 = l-2(bi+b 2 ) 

64 = t>2, 65 = Ol 



ai = 1/2 - ^3/28 
<J2 = 1/2 — ai 
0,3 = 0,2, Oi — a\ 



The 5-stage (8,4) method: Nig 



61 = 0.811862738544516 
02 = -0.677480399532169 

63 = 1/2 -(61 +62) 

04 = 03, 65 = &2j b 6 = °1 



a L = -0.00758691311877447 
a 2 = 0.317218277973169 

a 3 = 1 - 2(ai +a 2 ) 
04, = 02, 05 — a\ 



ods, LF, are compared with the second order NI±{&,2) 
methods. In the lower panels we compare the RKNqA 
methods against the (8,4) methods jointly with the best 
among the previous second order methods. 

For a weak nonlincarity, when the system can be con- 
sidered as a perturbed harmonic oscillator, we clearly ob- 
serve that the HO-split is superior to the F-split. In this 
case, with a relatively small number of Hermite functions, 
it is possible to approach accurately the solution, but 
this procedure has a limited accuracy which can deterio- 
rate along the time integration and depends on the initial 
conditions. In addition, the methods addressed to per- 
turbed problems show the best performance: The (8, 2)h 
method performs best among the compared when a rel- 
atively low accuracy is desired and the (8,4)# method 
takes its place for higher accuracies. 

Figure [5] shows the results for a = 1. It is qualita- 
tively similar to the previous case yet the HO-split does 
not outperform the plain F-split (jll.21|) as significantly 
as before. Nevertheless, it is important to observe that, 
again, the best result is obtained for the HO-split. Notice 
that a higher number of Hermite basis functions is neces- 
sary to achieve the same accuracy as the Hermite-Fourier 
decomposition. 

Figure [5] shows the results for a = 100. The HO-split 
cannot be expected to be particularly useful because the 
system is far from being a harmonic oscillator. From 
Fig. [2] we expect a great number of Hermite basis func- 
tions to be required for a sufficiently accurate expansion. 
The results in Fig. |6] demonstrate this rather intuitive 
expectation, i.e. almost negligible precision despite the 
large number of basis terms M = 150. Remarkably, the 
proposed HO decomposition does not show these limita- 
tions and reaches the precision of the F-split (jll.211) be- 
cause we are solving the harmonic potential exactly up to 
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FIG. 3. Error versus the number of exponentials in logarithmic scale for different splittings for the Leap-Frog method. 



spectral accuracy. For this problem, we observe that the 
(8, 2) method has the best performance when a relatively 
low accuracy is desired, the (8, 4) method shows the best 
performance for medium accuracies and the RKNqA is 
the method of choice for higher accuracies. 



B. Time-dependent harmonic oscillator perturbed 
by weak quartic anharmonicity 

We consider now a harmonic oscillator with time- 
dependent frequency and perturbed by a weak static 
quartic anharmonicity 



d I 



1 



p 2 + -uj 2 {t)x A )ip + 



1 4, 
SQ-X tp. 



(IV.2) 



We first consider the case co 2 (t) = A(l + ecos(wt)) with 
w = 1/2, A = 4, e = 0.1, eg = 0.01. As reference, we 
take a highly accurate numerical approximation as ex- 
act solution and restrict the spatial domain to [—20, 20] 
for all experiments in this subsection . We compare 
the Hermite-Fourier method with the plain Fourier split, 
since Hermite polynomials are not appropriate in a time- 
dependent setting. For fast oscillating systems and if high 
accuracy is needed, the two-exponential fourth-order ap- 
proximation of the harmonic oscillator (|II.13|) can be im- 
proved by taking, for example, a higher order Magnus 
expansion (|II.12|) . As we have seen, the solution of 



^' = Qp 2 +^w) 



u 



can be written as 
U(t,0) =. 



-i^(ux 2 +f3(xp+px)+jp 2 ) 



(IV.3) 



(IV.4) 



and we have considered, for example, a sixth-order Mag- 
nus integrator [17| to approximate the evolution operator 
for one fractional time step, a,ih, i.e. U(t, t+aih). This is 



equivalent to take in PV.4P t — aih and the parameters 
a, /?, 7 are given by: 

a = Yg(5wi + 8uj 2 + 5w 3 ) 

+ T§6-(x( w i + w l) + 8w 2 + ^1W 2 + u 2 lo 3 - ^WiW 3 ), 
P = a t hJ^(u} 3 - Wi)(i + %jg-(5wi + 8w 2 + 5w 3 )), 



7=1 



54 



(wi - 2lu 2 + uj 3 ) 



with u>i — u)(t n +Cih), i = 1, 2, 3 and c\ — 1/2 — v 15/10, 
c 2 = 1/2, C3 = 1/2 + \/l5/10, corresponding to a sixth- 
order Gaussian quadrature rule. The obtained operator 
is then decomposed according to Theorem 111.31 

The results are given in Fig. [7] and corroborate the 
superiority of the HO split. 

Another interesting example is given by an intense 
short pulse modeled via 



w(t) = Wq I 1 + 



At 



cosh 2 (B(t- 2)) 

Varying the parameters A and B, the pulse can be 
sharpened while keeping its time-average and hence its 
strength relative to the anharmonicity constant. Fig- 
ure [5] shows the results obtained for a relatively slow 
variation of the harmonic potential, for the parameters: 
wo = 4, A = 0.25, B = 2. Again, the advantageousness 
of the presented decomposition can be appreciated. It is 
already noticeable, that the error introduced by the time- 
dependence becomes dominant, and this effect increases 
for more rapidly varying potentials, e.g. for B> 1. In 
that case, higher order approximations of the Magnus 
expansion are necessary to maintain the benefits of the 
Hermite composition. 



V. CONCLUSIONS 

We have presented new Fourier methods for the nu- 
merical integration of perturbations to the time depen- 
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FIG. 5. Same as Fig. g]for a = l. 



dent harmonic oscillator which are useful for both both 
the Gross-Pitaevskii equation as well as for the linear 
Schrodinger equation. Fourier methods have shown a 
high performance in solving many different problems 
which can be split into the kinetic part and a remain- 
der that is diagonal in the coordinate space. We have 
extended the Fourier methods to perturbations of the 
time-dependent harmonic potential, and refer to them 
as Hermite-Fourier methods. They solve the linear 
Schrodinger equation with a time-dependent harmonic 
potential to the desired order using corresponding Mag- 
nus expansions and up to the accuracy given by the spa- 
tial discretization. These methods are fast to compute 
since FFTs can be applied and show a high accuracy 
when the problem is a small perturbation of the har- 



tend to perturbed harmonic potentials in linear quantum 
mechanics, c.f. section H*VB[ where it is straightforward 
to generalize the results to higher dimensions. 
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FIG. 8. Compare Fig. [7) The parameters used are wq — 
4, A — 0.25, B = 2 with a small anharmonicity eq = 0.01 



